Entanglement replication in driven-dissipative many body systems 
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We study the dissipative dynamics of two independent arrays of many-body systems, locally driven by a 
common entangled field. We show that in the steady state the entanglement of the driving field is reproduced 
in an arbitrarily large series of inter-array entangled pairs over all distances. Local nonclassical driving thus 
realizes a scale-free entanglement replication and long-distance entanglement distribution mechanism that has 
immediate bearing on the implementation of quantum communication networks. 



PACS numbers: 03.67.Bg, 42.50.Dv, 03.65.Yz, 42.50.-p 

Driving quantum systems to desired target states with very 
high fidelity is a central goal in quantum sciences and tech- 
nologies, in order to realize efficient and scalable devices be- 
yond the current state of proof-of -principle demonstrations. 
In pursuing this end, it has surfaced in recent years that the 
effects of noise and dissipation do not have to be necessar- 
ily detrimental in the realization of quantum coherent struc- 
tures (jJ-qII- The possibility of using suitably engineered ir- 
reversible dynamics to control quantum many-body systems 
has been discussed in a variety of settings, including driven- 
dissipative ultracold atoms in optical lattices [6], the asymp- 
totic realization of entangled states and quantum computa- 
tion in quantum spin models |7J,|8|], the dissipative control of 
trapped ions [9], and the steady-state entanglement of macro- 
scopic atomic ensembles [10]. On the other hand, ever since 
the formulation of the proposal for quantum repeaters 01 ll and 
the design of schemes for the implementation of remote quan- 
tum communication and distributed quantum gates [12], quan- 
tum networks have emerged as the strongest viable paradigm 
for the "quantum internet", i.e. the implementation of scal- 
able quantum computation and information processing sat- 
isfying the combined requirements of robustness, flexibility, 
multi-tasking and long -reach 11311 . A key ingredient of a quan- 
tum internet is the ability to hybridize, i.e. to interface het- 
erogeneous subsystems in a reliable and reproducible way. 
The strive towards the realization of such interfaces has been 
boosted by recent ground-breaking demonstrations of high- 
efficiency entanglement and state transfer between light and 
matter systems [fl4l — Tl6fl . and light-mediated teleportation be- 
tween remote nodes of a simple quantum network II 1711 . 



In this context, light-matter interfaces for the distribution of 
entanglement among network nodes, which exploit the robust- 
ness of irreversible dynamics, have been explored in several 



works Ill8l - l20ll . There it was shown that a reservoir of en- 
tangled light can drive distant matter systems into entangled 
states, thereby realizing an efficient transfer of entanglement 



from continuous- to discrete-variable systems. 

In the present work we show that when considering inde- 
pendent arrays of many-body quantum systems this mecha- 
nism amounts to the replication of the driving entanglement 
over many pairs of subsystems across the initially indepen- 
dent arrays. Specifically, we address the irreversible dynam- 
ics of two non-interacting chains of quantum systems simul- 
taneously driven, on one of their ends, by an entangled two- 
mode squeezed field (squeezed bath). The constituents in 
each array are coupled by nearest-neighbor linear interactions 
whose specific form is introduced below for different mod- 
els. The competition between the "entanglement pumping" 
process and the intra-array couplings results in a steady state 
consisting of a series of inter-array entangled pairs, each in- 
volving subsystems occupying corresponding sites in the re- 
spective chain [See Fig. UJ. Thereby, an arbitrary number of 
copies of identically entangled states is generated across the 
two arrays without violating fundamental constraints such as 
the no-cloning and the no-broadcasting theorems 02 ill . 

The replication mechanism works efficiently in different 
settings such as chains of harmonic oscillators or of spins. 
For pure harmonic resonators in the stationary state exactly 
N inter-chain pairs are formed that replicate the driving state 
independently of the size of the arrays. For two-level systems 
an ideal Einstein-Podolsky-Rosen (EPR) driving field creates 
exactly N Bell states across the two chains. 

To start, let us consider two chains of resonators, realiz- 
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FIG. 1 : A pair of independent arrays of linearly coupled quantum 
systems is locally driven by a two-mode entangled field. The ele- 
ments in each array are labeled by the indices je [1,N] (first chain) 
and je [N+ 1, 2N] (second chain). The steady-state inter-array entan- 
gled pairs are marked by dashed arrows. 
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ing two disjoint Jaynes-Cummings lattices J22 , 23 ], which can 
describe, in limiting cases, the physics of different condensed- 
matter systems ranging from spin chains, to boson or fermion 
lattice models. The two arrays are assumed equal (deviations 
from this condition are discussed below), and each consists 
of N single-mode cavities with equal resonance frequency 
and corresponding annihilation (creation) operators hj (St). 
Cavities belonging to the same array interact via nearest- 
neighbor linear coupling with strength rjj. Moreover, each 
cavity can interact resonantly with a two-level system (e.g. an 
atom in the cavity) with lowering (rising) operator <x ; (d"t), 
As illustrated in Fig. [TJ the elements of the first (second) ar- 
ray are labeled by indices j e [1,N] (j e [N + 1,2AT|). 
The two end cavities 1 and N + I are driven by a two- 
mode squeezed field. Including the dissipation of the cavity 
fields 124(1 , the master equation describing the system dynam- 
ics is p — —i[H c + H cs ,p] + Hup + -CsP- The unitary part 
of the evolution is ruled by the Hamiltonian H c +H cs with 
H c =Zy=i 1 rij{a\aj + i + a^ +j a N + j+ i + h.c.) describing the co- 
herent cavity dynamics and H cs = ~Z/j=\ gj(o~]®j + ^"If+Afr+j + 
h.c.) accounting for the interaction (with coupling gj) be- 
tween cavity j and its two-level system. The term X.d ac- 
counts for the dissipation of the cavities (at rate kj) and reads 
-CdP= 2/=i Kj{2djpa\-\a.aj,p}). Finally, £s accounts for the 
driving (at rate of the first-end pair of cavities ( 1 , N + 1 ) by 
the external two-mode squeezed field 



£,sp=2£ m(a xpa^+i +a n+ \pa i— SiOa/+ ip—pa i a^+i +h.c .) 
+ ^ £ [{h+Y){2a jpa\— {a\aj,p\)+h{2a\pcij— {Sj-at,p})] . 

;=l,iV+l 

The sum is over indices j-1 and j=N + 1 only, while h 
and in are related to the statistics of the driving two-mode 
entangled field: h is the same average photon number for 
both modes, m accounts for the inter-mode correlations, and 
m< y/h(h + 1), with equality holding in the squeezed vacuum. 
This effective model is based on the elimination of the de- 
grees of freedom of the reservoir (the driving field) in the 
limit of large squeezing bandwidth The entangle- 

ment in the driving field is the resource to be transferred via 
the replication mechanism. The state of the driving field is 

P s y=u inPT ul with (j in = e !^ifihl-KK) where K and hu 

are the field mode operators and p T a thermal state with n T 
average photons. The condition of large squeezing band- 
width corresponds to an almost constant squeezing param- 
eter, r(u>)~ro, over a sufficiently large range of frequencies 
around the cavity resonance. In this situation, the parameters 
characterizing the entangled driving field are n—nj+{2n.T + 
l)sinh 2 ro, ra=(«r+l/2) sinh(2ro). The entanglement is quan- 
tified by the logarithmic negativity £#= max[0, - log v_] with 
v-=2h+\-2m the smallest symplectic eigenvalue of the par- 
tially transposed covariance matrix for the two-mode field 
i25tl . The state is entangled iff v <1, which implies m>h. 

An exact analytical solution for the steady state is obtained 
if the arrays are driven by a two-mode squeezed vacuum (m - 
-\Jh(h + 1)), and X.d-0- To obtain the steady state in this situa- 



tion, we exploit the squeezing transformation U-®^ Uj^+j, 

with Uj^ + j=e ( ~ V> ' ro( - & ^ N *i~ & i &N *i ) , which maps the system into 
an equivalent one, whose density matrix p = U^pU satisfies 
the master equation p=-i[H c + H cs ,p] + -CsP = -Cp- The new 
dissipative term reads Xsp- lZj=i ,n+i <T [(2fljpaj-{flja>p})], 
and the transformed Hamiltonian for the cavity-atom 
interaction is H cs = ^J =1 gj[&^Cj(n) + al +j Dj(n) + h.c], 

with Cj(n) = V« + 1 &j + (-iy &if +j and Dj(n) = 

V« + 1 d~j+N + (-I) 7 This shows that, in the new rep- 

resentation, the arrays are in contact with a vacuum reser- 
voir and that each field mode interacts with two atoms at 
sites (j, N + j). It turns out that, regardless of the actual 
values of gj and r\j, V/' e \\,N\, the unique steady state 
is the pure state (which satisfies £\<p)(<p\ = 0) of the form 
\<p) = \4>) ® N j = \ |0> Q)j,N+j> i- e - the tensor product of the trans- 
formed modes' vacua with the atomic entangled state 



\4>y- 



i-4\i,i)j M+j +(-iy +1 c n \2,2) m+J ] 



(i) 



Here |1) and |2) indicate the ground and excited atomic states, 
and c„ - V n /(2« + 1)- Due to the destructive interfer- 
ence between transitions amplitudes involving the atomic pair 
(J, N + j) that is coupled to the same mode, state \<p) is such 
that the atoms are decoupled from the field. Moreover, it is 
not affected by dissipation because the field modes are in their 
vacuum state. Therefore, during the dynamics, population ac- 
cumulates, eventually pumping the system into the entangled 
state Eq. ([TJ). Going back to the original representation (by in- 
verting the transformation U) also the field modes become en- 
tangled in inter-array two-mode squeezed vacua for each pair 
(j,N + j): Ujjy+j\Q,0)jjf + j. All inter-array field-pairs have the 
same entanglement of the input driving field, thus realizing 
a perfect entanglement replication mechanism. On the other 
hand, the entanglement of all inter- array atomic pairs is the 
same as that discussed in Refs. [ 18-20J] for a single atomic 
pair, but with the essential difference that it is now exactly 
replicated across all the pairs. This is the main result of this 
Letter, from an ideal, infinitely entangled state of the driving 
field one obtains by engineered dissipation an arbitrary num- 
ber of EPR field pairs and Bell states of the atomic pairs. In 
general, the entanglement of the pairs is limited only by the 
amount of entanglement of the driving field. Moreover, as 
will be shown below, this result is rather general as it holds 
valid also for spin chains and arrays of harmonic oscillators. 

We will now study the effects of a non-negligible thermal 
nature of the driving field, and of other sources of dissipa- 
tion and noise. We consider first the limit in which the model 
reduces to two chains of harmonic oscillators, i.e. when the 
atoms are not present {gj — V j). In this case an exact analyt- 
ical solution is found also if the external field is not perfectly 
squeezed, m< y/h(h + 1). We still assume that = 0, and 
we find that, in the squeezed representation, the steady state 
of each cavity is thermal, p~ , with mean occupation number 
hj. In the anti-transformed representation, this corresponds 
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FIG. 2: £^ DV ' as a function of (a) the pair-site label (j, N + j) (with N=20 and k^ko V j), (b) /V (with Kj=0. It], V;'), (c) « (with m= + 1)), 
(d) m (with n=l), and (e) k n =k 2 n (with a- j¥jV ,2jv=0)- In all cases rjj = 77, Vj, and ^=77. The remaining parameters in (a), (b) and (e) are «=1 
and ih= ^lh{h + 1); in (c), and (d) they are Kj=0.1t] and N=10. The insets indicate the pair (J, N + j) corresponding to each line. In (e) the 
dash-dotted curve corresponds to all pairs (j, N+ j) for j e [2, N-l]; These results are independent of N and have been verified numerically for 
arrays of size up to N = 30. The solid thick (gray) lines report the entanglement of the driving field which is equal to the entanglement of each 
pair when k ; = 0, V j. 



to a two-mode squeezed thermal state for each pair of field 
modes (j,N + j) that reads: U jj+nP^ ® Pt +N)u Ij + n- The 
corresponding steady-state entanglement is the same as that 
of the driving field, regardless of j, N, h and m. Therefore, 
the exact replication of the driving field entanglement takes 
place also in this case. When the other sources of dissipation 
described by are included, the steady state of the system 
can be determined numerically, and the logarithmic negativ- 
ity &n[], k] of any pair (j, k) of cavity fields is obtained from 
the corresponding covariance matrix ll25ll . Quantitatively, we 
study the logarithmic negativity normalized to unity, defined 
as E { ™ v) [j, k]=S N [j, k]/(l + S N [j, k]). 

Most of the results to follow are obtained for a reservoir 
with h = 1, such that the corresponding entanglement is rela- 
tively small. Remarkably, even in this strongly non-ideal sit- 
uation, the replication mechanism is significantly resilient to 
the added noise. As shown in Fig. [2] (a), the entanglement 
decreases with the decay rate of the cavities. At fixed decay 
rate, the largest E^ m ' is achieved by the pair {\,N + 1 ) that is 
directly coupled to the driving field. The entanglement of the 
other pairs decreases moderately with the distance from the 
driven pair and exhibits a weak revival for a few pairs at the 
opposite end of the arrays. Fig. |2](b) illustrates how the entan- 
glement mildly decays with the size of the arrays, remaining 
nonvanishing up to large values of N. Hence, the entangle- 
ment replication mechanism exhibits a notable robustness in 
the presence of losses. The dependence of the entanglement 
on the statistics of the input field is shown in Fig. [2] (c) and 
(d). When the driving is a squeezed vacuum, its entanglement 
increases with h [gray line in Fig. [2] (c)] and reaches unity 
asymptotically as «— »oo. For lossy cavities, the entanglement 
saturates to a value smaller than unity that depends on the 
pair being considered. The entanglement distributed through 
a squeezed thermal state is reported in Fig. |2](d) showing that 
£ , ^" v) is nonvanishing for all values of m for which the driving 
field is entangled (m > h). When only the end cavities are 
open {kj^.in = 0), the pairwise entanglement is minimum 
at km = i<2N-r] for all pairs (j, N + j) except for pair (Af, 2N) 
whose entanglement instead decreases monotonically with 
[See Fig. |2](e)]. As increases the coherent coupling be- 



tween the last cavity of each array and the neighboring one is 
progressively inhibited. At large values of each of them is 
effectively decoupled from the rest of the system, whose en- 
tanglement is thus restored to the value of the non-dissipative 
case. Moreover, the field leaking out of the last pair of cav- 
ities is entangled as well [27] and even equal to that of the 
driving field for some frequencies [27]. This feature allows 
for the re-usability of the transferred entanglement for net- 
working protocols. So far we have discussed results obtained 
with homogeneous couplings rjj=i]. Analogous results hold 
even with intra-array patterns of inhomogeneous couplings, 
as long as the two arrays remain equal. Asymmetries between 
the arrays reduce the inter-array entanglement, but the repli- 
cation mechanism remains valid as long as they are not too 
strong. This is shown in Fig. [3] (a), obtained for random cou- 
plings rjj = t]q + gj, with j e [l,2N], where are zero-mean 
random variables uniformly distributed in a range A£. 

When each cavity interacts with a two-level atom we can 
study the entanglement properties of the atoms by approxi- 
mating the system with an effective spin model. We focus 
on the weak coupling limit, such that the couplings gj be- 
tween the atoms and the cavities are sufficiently small |EtTi 
and we can adiabatically eliminate the cavity fields to find a 
closed equation for the atoms. The resulting spin model ex- 
hibits non-trivial long-range interactions and collective decay 
of the spins, as reported in detail in the Supplementary Ma- 
terial B27I1 . Here we discuss the results relevant for the cor- 
responding steady state. Let us consider the logarithmic neg- 
ativity [j, k] = log 2 llp^Jlli of the state pjk of the atomic 
pair (j, k), where || ■ ||i is the trace norm and PT stands for par- 
tial transposition. The entanglement properties of the atoms 
are similar to those of the free cavity fields. However, at vari- 
ance with the latter case, '[j, k] is sensitive to the statistics 
of the driving entangled field and decreases more rapidly with 
decreasing m as illustrated in Fig. 0(b). 

The effective spin model with long-range interactions 
can be compared with the case in which two indepen- 
dent spin chains with XX short-range interactions are cou- 
pled on one end to the driving field. As shown in Fig. [3] 
(b) and (c), one obtains very similar results. The master 
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FIG. 3: (a) for a model with random couplings as specified in 
the text. The curves are obtained averaging the result over 500 real- 
izations. For each value of A£ the vertical bars represent the interval 
between the realizations of maximum and minimum entanglement. 
The other parameters are N = 10, h = 1, m = sjh^i + 1), f = rj , 
and Kj = 0.02;7oV j. (b)-(c) Comparison between the logarithmic neg- 
ativity for atoms in cavity arrays, Eft (b), and for spins in XX spin 
chains, (c), as functions of m for h = 1. The remaining param- 

eters are kj=0 Vy, N=3, (=tj, and g=O.Ol77 for the atoms, and N = 3, 
and Jj = y V/ for the spins. The insets specify the correspondence 
between curves and pairs (j, j + AO- 



equation for this case reads p = -i[H s ,p] + -Csp, with 
« s =lZyr i 1 Z*=o,iv-/j(^ + ^ +i+1 +^ + ^ +i+1 ), where Jj is 
the spin-spin coupling, and &*' y are the Pauli spin operators. 
The effect of the driving field is described by 



-CsP 



2m(&ip& N +i+fr N+ ip&i-&i& N +ip-p&i& N+ i+h.c) 



+ [(n+l)(2o-7P^-{o-]>j,p})+n(2o-]po- ; -{o- j a-],p})] , 

j=l,N+l 

with cr, (cr[) the spin lowering (rising) operator. While in 

J j 

the cavity-atom system the effective spin-spin interactions are 
long range [27], here we deal only with local ones. Neverthe- 
less, entanglement replication continues to hold. Indeed, the 
stationary state of the system for m= ^h{h + 1) can be evalu- 
ated analytically and coincides with that of Eq. (Q3, where |1) 
and |2) now denote, respectively, the spin up and spin down 
states. Finally, we observe that the similarity of the steady- 
state entanglement properties in the two systems holds even 
when the driving field has a nonvanishing thermal component, 
as shown in Fig. 0(b) and (c). This result shows the generality 
of the entanglement replication mechanism which is largely 
independent of the specific physical realization. 

In conclusion, we have discussed a scheme realizing the 
replication of entanglement, based on the interface of a driving 
two-mode entangled field with two distant and independent 
dissipative many-body systems. The replication mechanism 
works efficiently both for arrays of discrete- and continuous- 
variable systems. Since the phenomenon occurs in the steady 
state of the irreversible driven-dissipative dynamics, it ex- 
hibits an intrinsic robustness against the detrimental effects 
of noise. We have highlighted the roles played by quantum 
interference and the competition between dissipation, driving, 
and interactions in producing such a steady state. The corre- 
sponding entanglement is robust against deviations from ideal 
conditions including a nonvanishing thermal component of the 



driving field, asymmetries between the arrays, and decay of 
the cavity fields. Ideally, the replication mechanism yields an 
arbitrary number of maximally entangled pairs and is scale- 
free in the sense that it is independent of the actual length of 
the arrays. Thus, it is a potentially valuable resource for re- 
mote quantum communication and distributed quantum com- 
putation 11121 1 1311 that could be combined with other driven- 
dissipative strategies for the realization of scalable quantum 
networks B28I1 . Seen from a different viewpoint, this scheme 
implements a protocol of long-distance entanglement distribu- 
tion (29, 30j] and nested entangled-pair production ll3~lll . two 
key tasks for quantum networking, achieved via the interac- 
tions intrinsic in many-body systems. 

The outlined scheme is general and flexible enough to 
find application in many systems which effectively realize 
chains of harmonic oscillators or spins, such as cavity/circuit- 
QED I32L 13311 . arrays of optomechanical systems, trapped 
ions, or ultracold atoms in optical lattices. The mechanism 
could be verified with arrays of coupled resonators, recently 
produced in photonic crystals ll34ll35ll . which realize chains of 



linearly coupled harmonic oscillators. In Ref. [34] the cavities 
are almost resonant and they interact with nearest-neighbor 
couplings of strength within the range ~ 60 - 2000 GHz. 
These values can be tailored by selecting the distance between 
the cavities. The reported cavity line-width is of the order of 
~1 GHz. These parameters are consistent with those discussed 
in our analysis. However, the broadest squeezing at the wave- 
length of the resonators of Ref. fl34ll (~ 1 .5 pm) has a band- 
width of about ~2 GHz 13611 . This value is still relatively small 
and does not well satisfy the broadband condition assumed 
throughout our work. Nevertheless, larger squeezing band- 
widths and photonic-crystal nano-cavities with weaker decay 
rates are expected to be realizable in the near future 113 6i 13711 . 
thus matching the required condition. On the other hand, the 
currently available experimental situation might already suf- 
fice for testing the entanglement replication mechanism. In- 
deed, a relevant theoretical question, which deserves further 
investigation, is whether entanglement replication holds also 
for driving squeezed fields of finite bandwidth. 
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SUPPLEMENTARY MATERIAL 



LOGARITHMIC NEGATIVITY OF THE OUTPUT FIELD 

If the last pair of cavities of the two arrays are open, then 
the field leaking out of the last pair of cavities is entangled. 
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The corresponding logarithmic negativity En ' [See Fig. [4] 
(a)] is maximum at intermediate values of kn- At small values 
of kn the number of photons leaking out of the cavities is too 
low, and so is the associated entanglement. Similarly, at large 
k n dissipation is too strong for the build-up of entanglement 
in the output fields (See below for the detailed evaluation of 
E N ° ). Remarkably, in this situation, E N ° can reach val- 
ues very close to those of the driving field, thus demonstrat- 
ing the effectiveness of the scheme and the re-usability of the 
transferred entanglement for networking protocols. 

The cavities in the arrays can emit photons into the con- 
tinuum of modes of the electromagnetic field of the environ- 
ment. Therefore, the output field is made of a continuum of 
frequencies. In order to determine the logarithmic negativity, 
E^^ioS), of the frequency components of the output field, we 
have to evaluate the spectrum of the covariance matrix. There- 
after, the value of E ( ° ut \ai) is obtained by applying the defini- 
tion of the logarithmic negativity to the spectral components 
of the covariance matrix. 

An example of E^"'\aj), corresponding to the parameters 
for which the entanglement of the output field reaches, for 
some frequencies, a value very close to that of the driving 
field, is shown in Fig. H](b). Here maxima of the entangle- 
ment are found in correspondence of the frequencies of the 
normal modes of the arrays. The figure Fig. |4](a) illustrates 
the behavior of En as a function of kn, evaluated in terms of 
the maximum value of the spectrum (of), for each value 



of the decay rate kn, that is En = max 

In general, the covariance matrix of the output field can be 
expressed as 

r (ou,) (w)=i [0yt (o, " ) (w)0 T + &m (out \co) T & T ] . (2) 

where the elements of the 4Nx4N matrix are 0^ = 5/,2/c-i + 

<5;,2(t-4A<-i + i (f>j,2k _ <5/,2yt-4iv), an d J?T (o " f) (w) is the spectrum of 
the correlation matrix of the output field operators defined as 



with 



CO 



(3) 



(4) 



where a,B e {+, -}, and we use the definitions a + = a . , 

1 ' J ' j out join 

and aj gul = a ]out for the creation and annihilation operators of 
the output field Q. 

The correlation functions of the output field, in Eq. 
can be evaluated by means of the input-output formalism [38] 
which allows to express the correlation functions of the output 
operators in terms of the correlation functions of the system 
operators. Exploiting this formalism, one finds 



J 



rt ol "\oS) = - 



' K [(M" + \uyTf x k~ + A - r (M- - iwl)- 1 ] K K [(M" + iwl)- 1 A+- r + AJ" 



(M + 



j1) _1 ]K- 



K [(M + + iwl)- 1 A+- + AJ-(M~ - iwl)- 1 ] K K [(M + + ia»l)" 1 A+ +r + AJ + (M + - iwl)" 1 ] K 

I 



where K is the 2N x 2N diagonal matrix with elements Kjj = 
-JFi, the matrices M" are the matrices of the coefficients in the 



system of equations for the evolution of the averages of the 
cavity field operators j,{a a j) = Yik M^S"), and the elements 

of the matrices A^f are the steady-state correlation functions 
of the cavity field operators, defined as 



(5) 




10 10 10 10 10 

K H 



The elements of the matrices M a are easily evaluated, whereas 
the matrices A^ 8 can be computed numerically solving the set 
of equations for the correlation functions whose form is found 
using the master equation for the system dynamics in the main 
text of the present work. 



FIG. 4: (a) E N as a function of k n =k 2 n (in units of 77). (b) E ° as 
a function of a> (in units of 77) for k n = kin = OArj. The frequency oj 
is relative to the resonance of the cavities. The vertical dotted lines 
indicate the frequency of the normal modes of the arrays. In both 
panels the remaining parameters are Kj^ Nt 2N-0, N=10, £=0.5rj, h=l, 
m = \2. The solid (gray) lines indicate the entanglement of the 
input squeezed vacuum. 
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EFFECTIVE SPIN MODEL FOR THE TWO-LEVEL 
SYSTEMS DYNAMICS 

Here we consider the model described in the main text of 
the present work, with homogeneous couplings rjj = tj and 
gj = g for all j. When the time scale for the two-level 
system dynamics (T at ) is much larger then the time scale 
for the cavity-fields dynamics (T cay ), then we can adiabat- 
ically eliminate the cavity fields, thereby obtaining an ef- 
fective spin model for the dynamics of the two-level atoms. 
The time scale for the atoms dynamics can be estimated as 
~ l/g V(") + 1, where («) stands for the average cavity- 
photon number, whereas the time scales for the fields dynam- 
ics is determined by the eigenvalues of the matrix M" 
defined in Sec. „ In particular the time scale for the field dy- 
namics is set by the smaller eigenvalue T c - dv ~ 1 /min 

Henceforth, the master equation for the system dynamics 
can be rewritten as 



p = £oP + t\p, 



(6) 



where £ p- - i[H c ,p]+£ s p+£ D p and £ x p- - i[H cs ,p]. At 
lowest order in the atom-field coupling strength, the master 
equation p s = Tr/,-^ [p] describing the dynamics of the two- 
level atoms only takes the form 

dt -Yr field (Xie^'Xip^ ® Ps ) , (7) 

where p{' t e,d is the steady-state of the field in absence of the 
interaction with the atoms. This expression can be recast as 

AN 

ps- ^ [o-j (Tkj + Tj,k)p s &k-&jTj,kd-kp s -p s d-jTkjd- k ] (8) 

j,k=\ 

where &j = &'. for j < 2N and &j = Q-j-in otherwise. We 
have introduced the AN x AN matrices T and T with elements 



AtTx fidd {a f £a "Rk) (0 = T,T) 



(9) 



with <R k = akPst {Rk = p sl a k ) for O = T (O = T). Here dj = 
a,j for j < 2N and dj = d\_ 2N otherwise. Both T and T can 

be expressed in term of the matrices M ff and A? 8 , introduced 
above in Section n of this supplementary material, as 



T = 



T 



(M-)- l A, 



v (M + )- 1 A, 



(M-) 
o ( M+ ) 



-1 A ++ 



A + (M + )- 



a ( ; + (M + )- 



(10) 



Equation ([8]l and the matrices in Eq. ( TTOb have been used for 
the numerical evaluations presented in the discussion of the 
atom-cavity model. Eq. (O describes a non-trivial spin system 
where both the spin-spin coherent interactions and the dissipa- 
tion mechanism can be long-range. An example where such 
effective spin model can be studied analytically is found for 
-Co = 0, as seen in the next Subsection. 



Effective spin model for £ D = 

When £d = the effective master equation takes the form 

AN 



Ps = 7f i 2 (T j^KjP'' a ' k ~ a J^fk^kPs ~ PsCT j^O-k] 

j,k=\ 
AN 



j,k=i 

(id 

where £ e {even, odd) distinguish between the case in 
which N is even or odd, the parameters are 



J — i yeven — K , , y odd ~ > 

77 rf k 



(12) 



and the AN x AN matrices of coefficients ?(® and J/®, can be 
expressed as block matrices 



(l+«)X f ^ 



-nXf 



— «X^r 

mW ? (1 + «)Y f 



(13) 



mW* 



Here the missing blocks are null matrices, and are NxN 
matrices whose elements are 

N 

0^eren)j,k = ^]( _ 1)" +1 \d j,2m8 j,k+2n-l +^k,2m^j+2n-l,k\ . 
N 

(X dd)j,k =^]( _ 1)" +1 [6j,2m+l6j,k+2n-l+6k,2m+l5j+2n-l,k\ , 



n,m= 1 



(Yeven) j,k~ 2u j,2m^ jM2n + <>k,2m<) j+2n,k\ 

n,m=l 
N 

m-1 

N 

(Yodd)j,k= X (-1)" \dj,2m-\6j,k+2n + 5k,2m-\5 j+2n,k\ 
n,m= 1 
N 



and W f = (y ? + i^Xfjz with Z M = (-1)' ; (V U . 

The first term in Eq. (fTTT l describes the coherent interaction 
between the spins, while the second one accounts for the dis- 
sipation. The coherent part does not couple spins belonging 
to different arrays, and the spins in each array are coupled ac- 
cording to the structure defined by the matrix in Eqs. ( TBI : 
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the indices of the nonvanishing entries in these matrices corre- 
spond to the indices of the coupled spins. The incoherent part, 
on the other hand, couples both spins from the same array and 
from different arrays, according to the pattern defined by the 
matrix J/®, in Eq. (fT3T >. 
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